The prior art for atmospheric correction of spectral imagery (that is, the retrieval of surface reflectance spectra from measured radiance spectra) at wavelengths where the source of light is the sun is embodied in various methods described in the literature. The simplest method is the “Empirical Line Method”. It assumes that the radiance image contains some pixels of known reflectance, and also that the radiance and reflectance values for each wavelength channel of the sensor are linearly related; therefore, the image can be converted to reflectance by applying a simple gain and offset derived from the known pixels. This method is however not generally applicable. In-scene known reflectances are often not available, and the linearity assumption, which presumes a uniform atmosphere throughout the scene, may not be accurate. Therefore, an alternative method based on a first-principles radiation transport model is required. The phenomenology and methods differ somewhat for images of land and water; the focus here is on methods appropriate for scenes consisting mainly of land.
For hyperspectral imagery, which contains many tens to hundreds of wavelength channels, the most commonly used atmospheric correction computer code at the present time is ATREM [Gao et al., 1996]. The inputs to ATREM consist basically of the center wavelengths and widths of the sensor channels, the location and time of the measurement(which defines the sun angle), the sensor altitude and surface elevation, a specified aerosol/haze type and amount (as defined by a visible range), and a specified default model atmosphere (that is, one of a choice of embedded atmospheric descriptions, such as a seasonal/latitude model or the US Standard model, etc.). The sensor view is assumed to be in the nadir direction. Another atmospheric correction code for hyperspectral imagery has been developed by Green et al. [1996] for processing data from the AVIRIS sensor. Green's method requires additional input from a calibration image, and is therefore not applicable to general hyperspectral data. Recently, a method has been developed by Spectral Sciences, Inc. (SSI) and the Air Force Research Laboratory (AFRL) [Adler-Golden et al., 1999] that combines the essential features and inputs of ATREM with a more advanced treatment of the radiation transport physics similar to that found in the multispectral atmospheric correction code ATCOR2 [Richter, 1996]. The radiation transport model used in the SSI/AFRL method is the most recent version of the MODTRAN band model-based code, currently MODTRAN4 [Berk et al., 1998]. ATCOR2 and Green's method are also based on versions of MODTRAN. This method is suitable for off-nadir as well as nadir viewing. In contrast, ATREM uses the 6S code [Vermote et al., 1994] to model atmospheric scattering and accounts for molecular absorption using a transmittance factor derived from a separate model.
For technical details of the general first-principles atmospheric correction methodology, one may start with a standard equation for spectral radiance L* at a sensor, as associated with a given surface pixel,L*=Aρ/(1−ρeS)+Bρe/(1−ρeS)+L*a  equation (1)which applies for UV through short wave infrared (SWIR) wavelengths, where thermal emission is negligible. Here ρ is the pixel surface reflectance, ρe is an average surface reflectance for the pixel and the surrounding region, S is the spherical albedo of the atmosphere (i.e., the atmospheric reflectance for upwelling radiation), L*a is the radiance backscattered by the atmosphere, and A and B are coefficients that depend on atmospheric and geometric conditions. The first term in Equation (1) corresponds to the radiance from the surface that travels directly into the sensor, while the second term corresponds to the radiance from the surface that is scattered by the atmosphere into the sensor. The distinction between ρ and ρe accounts for the “adjacency effect” (spatial mixing of radiance among nearby pixels) caused by atmospheric scattering. The adjacency effect correction may be ignored by setting ρe=ρ.
For a specified model atmosphere the values of the parameters A, B, S and L*a in Equation (1) are determined from the radiation transport model. For example, in the SSI/AFRL method these values are derived empirically from the standard MODTRAN outputs of total and direct-from-the-ground spectral radiances computed at three different surface reflectance values, 0, 0.5 and 1. In ATREM the parameters are derived directly from the 6S code and the molecular absorption model. In both methods the viewing and solar angles of the measurement and nominal values for the surface elevation, aerosol type and visible range for the scene are used as inputs.
The water vapor profile is not generally known for the scene in question, and furthermore it may vary across the scene. To account for this, the calculations to determine A, B, S and L*a are looped over a series of varying water column amounts, then selected wavelength channels of the image are analyzed to retrieve an estimated column water vapor amount for each pixel. For example, the SSI/AFRL method gathers radiance averages for two sets of channels, an “absorption” set centered at a water band (typically the 1.13 micron band) and a “reference” set of channels taken from the edge of the band. A 2-dimensional look-up table (LUT) for retrieving the water vapor from these radiances is constructed. One dimension of the table is the “reference” to “absorption” ratio and the other is the “reference” radiance. The second dimension accounts for a reflectance-dependent variation in the ratio arising from the different amounts of absorption in the atmospherically-scattered and surface-reflected components of the radiance. In ATREM this ratio variation is not taken into account, so a 1-dimensional LUT based on the ratio alone is used.
After the water retrieval is performed, Equation (1) is solved for the pixel surface reflectances in all of the sensor channels. The SSI/AFRL and ATCOR2 codes use a method described by Vermote et al. [1997] in which a spatially averaged radiance image L*e is used to generate a good approximation to the spatially averaged reflectance ρe via the approximate equationL*e=(A+B)ρe/(1−ρeS)+L*a  equation (2)The spatial averaging is performed using a point-spread function that describes the relative contributions to the pixel radiance from points on the ground at different distances from the direct line of sight. The SSI/AFRL code approximates this function as a wavelength-independent radial exponential, while ATCOR2 uses a simple square array of pixels. ATREM ignores spatial averaging entirely and sets ρe=ρ. This amounts to neglect of the adjacency effect, and leads to considerable errors in scenes containing substantial amounts of aerosols s(including haze) and/or strong brightness contrasts among the surface materials.
A procedure analogous to the water vapor determination can be used to retrieve an approximate scene elevation map from a hyperspectral image. Here the radiation transport calculations are looped over elevation rather than water vapor concentrations and an absorption band of a uniformly mixed gas such as O2 or CO2 is interrogated. Reasonable relative elevations have been obtained from AVIRIS images using the 762 nm O2 band (as a default) or the 2.06 μm CO2 band.
In the prior discussion it has been assumed that the quantity of aerosol or haze in the scene to be analyzed has been adequately estimated. Both ATCOR2 and the SSI/AFRL code incorporate methods for retrieving an estimated aerosol/haze visible range from one or more “reference” surfaces in the scene where the reflectance is known for some wavelength channel(s). Best results are obtained using short (i.e., visible) wavelengths and either a very dark surface, such as green vegetation or deep calm water, or a very bright surface, such as a white calibration target that is large enough to fill at least one pixel. In the SSI/AFRL method the aerosol/haze retrieval operates as follows. The calculations to determine A, B, S and L*a are carried out over a bandpass containing the selected sensor channel(s), but instead of looping over a series of water vapor values as described previously they are looped over a series of visible range values, e.g. 200, 100, 50, 33, 25, 20 and 17 km, that are evenly spaced in their reciprocals (optical depths). This process is undertaken only over spectral ranges or bandpasses where water vapor absorption has little impact. The user chooses the reference pixels and assigns them each a mean reflectance value for the selected channels. A 2-dimensional LUT relating the visible range reciprocal to the pixel radiance and the spatially averaged radiance is then constructed. The visible range is reported for each reference pixel by interpolating from the LUT, and if desired the results can be combined to form an average for the scene. It should be noted that the above method, unlike many others, takes the adjacency effect into account. When the reference pixels are taken from small areas in an image, such as calibration panels or isolated patches of vegetation, the adjacency effect correction becomes critical for an accurate result.
Limitations of the Prior Art
The prior art methods for atmospheric correction of spectral imagery over land work adequately for many but not all conditions. In particular, they assume cloud-free conditions, the presence of at least one material in the scene with a known reflectance at a visible wavelength, an accurate and consistent wavelength calibration throughout the image, and sufficient computing time to perform tens of mathematical operations per image pixel per wavelength channel. As described below, problems arise when these conditions are not met.
Cloudy conditions. Clouds pose several problems for atmospheric correction. Not surprisingly, pixels containing either a cloud or a cloud shadow are not converted to true surface reflectance by the prior art methods. Less obvious is the fact that clouds can adversely affect the reflectance retrieval accuracy even in cloud- and shadow-free parts of the scene. This is because the clouds contaminate the spatially averaged radiance L*e used in Equation (2) to generate the spatially averaged reflectance ρe. According to theory, ρe should account only for reflecting material that is below the scattering atmosphere. Clouds can be embedded throughout the molecular scattering layers, anywhere within the troposphere, 0-12 km. which are important sources of radiance at the very shortest (blue-violet) wavelengths. However, clouds usually lie above the aerosol and haze layers (approximately 0-2 km for the boundary) that dominate the scattering in the rest of the spectrum. Since clouds are typically brighter than the terrain, they can lead to an overestimation of ρe, and hence underestimated surface reflectance retrievals. Therefore it is desirable to identify and “remove” the clouds prior to the calculation of ρe.
Absence of accurately known surfaces. For the purpose of aerosol/haze amount retrieval, vegetation, water, or other dark surfaces can frequently be identified in a scene. However, reflectance values for these surfaces at appropriate wavelengths are often not known to within the accuracy needed (around 0.01 reflectance units or better) for an accurate retrieval. Even with “calibrated” surfaces the reflectance may not be known to within this accuracy because of measurement complications caused by non-Lambertian bidirectional reflectance distribution functions. Therefore, an aerosol/haze retrieval method is needed that is less prone to reflectance estimation errors.
Wavelength calibration error and non-uniformity. The prior art methods presume that the center wavelengths and shapes of the channel response functions (“slit” functions) are correctly specified by the input values. In practice, the actual wavelength centers in hyperspectral sensors can differ from their presumed values by a sizable fraction of the channel width. This leads to large errors (up to tens of percent) in the retrieved reflectance spectra near atmospheric absorption features. Prior art methods do not account for these wavelength errors except in so far as the user may attempt to model them by trial and error, which is a tedious and time-consuming process. In addition, in sensors with a two-dimensional focal plane array there is often a significant non-uniformity of the wavelength calibration across the image; this may be due to curvature of the image of the spectrometer slit on the focal plane (referred to as spectral “smile”), misalignment of the focal plane and spectrograph, or other effects. Prior art atmospheric correction codes do not account for these sources of wavelength non-uniformity.
Computing time constraints. For a typical hyperspectral image containing several hundred or more spectral channels and hundreds of thousands of pixels or more, the speed of the atmospheric correction is fundamentally limited by the mathematical operations required to generate the reflectance values for each pixel and channel from the Equation (1) parameters. The prior art methods (e.g, ATREM, SSI/AFRL), which account for pixel-specific values of water vapor and possibly other quantities such as ρe, require tens of operations per pixel-channel. Many of the operations are consumed in interpolating to find the appropriate A, B, S and L*a parameters for each pixel, and even more would be required to account for spectral “smile”, for example. A more efficient procedure is needed for high-speed (near-real-time) atmospheric correction.